MODULE COMFUNCS
!This module contains common functions 
!such as human capital prduction
!period utility function  

  USE GLOBVAR
  USE NRUTIL, only : NR_BIG2, NR_SMALL, NR_SMALL2, ARRAY2D, NR_EPS
  USE INTERPOL, only : locate, ARRAY1D_INIT  !, linear_inte2 
  use INTEGRATION, only : Gauss_Hermite
  use MATRIX, only : Cholesky, Matrix_Inverse

  IMPLICIT NONE
CONTAINS
  
  !assign parameter values from the vector. called in all nodes  
  subroutine sub_parse_parameters()
      implicit none
      integer(4) :: i, j, ii
      real(8) :: rtemp, rtemp2

      pa_epsilon = gr_Vpar(1)
      
      if (gcI_robust == 21) then !no H depreciation at all
          psigma = 0.0d0        !H depreciation
      else
          psigma = gr_Vpar(2)     !H depreciation
      endif
      
      if (gcI_robust==12 .OR. gcI_robust==13) then  !exo or lbd Hone
          palpha_I = 0.0d0
          palpha_H = gr_Vpar(4)
      else
          palpha_I = gr_Vpar(3)
          palpha_H = gr_Vpar(4)
      endif
      i = 5
      pXI_sd = gr_Vpar(i)  !5
      i = i + 1  !6
 
      !==============================================================
      ! assign the heterogeneity to proper group (defined by color_mpi)
      ! Heterogeneity in pi (pi*(I^alpha)*(H^beta))  and a0
      ! Cholesky decomposition (transformation)
      ! X=mean_x+sigma_x * Z1
      ! Y=mean_y+sigma_y * (rho*Z1+sqrt(1-rho^2)*Z2)
      !==============================================================
      if ((N_het_a0==1) .and. (N_het_pi==1)) then
         pa0 = gr_Vpar(i)    !6
         ppi = gr_Vpar(i+1)  !7
         pinitHmean = gr_Vpar(i+2)  !8
         pinitHsd = gr_Vpar(i+3)    !9
         i = i + 4           !10
      else
         pa0 = gr_Vpar(i) + gr_Vpar(i+2)*p1_Vhet_x(gI_het_group_p1(color_mpi),1)
         ppi = gr_Vpar(i+1) + gr_Vpar(i+3)*(gr_Vpar(i+4)*p1_Vhet_x(gI_het_group_p1(color_mpi),1) &
                 & + sqrt(1.0d0 - gr_Vpar(i+4)**2)*p2_Vhet_x(gI_het_group_p2(color_mpi),1))
         i = i + 5      !11
         pinitHmean = gr_Vpar(i)
         pinitHsd = gr_Vpar(i+1)
         i = i + 2     !13
         if (N_het_H0>1) then
            pinitH_a0 = gr_Vpar(i)   !13
            pinitH_pi = gr_Vpar(i+1)   !14
            i = i + 2  !15
         else
            pinitH_a0 = 0.0d0
            pinitH_pi = 0.0d0
         endif
      endif
      ppi = max(0.001d0, ppi)
      
      if (gcI_CF == 4) then
         pa0 = pa0 + 0.5d0*(pa_epsilon**2)
         pa_epsilon = 0.0d0
      endif
      
      pb1 = gr_Vpar(i)  !(15)
      i = i + 1    !i=16 now
      
      if (N_pbequest == 2) then
          pb2 = gr_Vpar(i)  !(16)
          i = i + 1    !i=17 now
      endif
      
      if (gcL_searchBC) then
          gr_BC_Ratio = gr_Vpar(i)   !16 the ratio of borrowing constraint (0,1)
          sgridAnrlb%a = gr_BC_Ratio * sgridAnrlb_BC%a
          do j = N_0 + 1, N_T
              ii = j - N_0 + 1
              sgridA%a(ii,1:VN_A(j)) = sgridA_BC%a(ii,1:VN_A(j)) + sgridAnrlb%a(ii)
          enddo !j
      else
          petac = gr_Vpar(i)  !(16) petac !consumption, CRRA, 1.0d0 - (-3.0d0)
      endif
      
      pphi(1:N_phi_t+1) = gr_Vpar(i+1:i+N_phi_t+1)  !consumption utility shifter (17-20)
      
      pa1 = gr_Vpar(i+N_phi_t+2)     !taste of leisure when spouse not working (21)
      pa2 = gr_Vpar(i+N_phi_t+3)     !taste of leisure when spouse working   (22)
      
      i = i + N_phi_t + 4   !23
      if (IC_L_nonseparable >= 1) then
         pphi(N_phi_t+2) = gr_Vpar(i)  !pphi5: C_L nonseparable
         i = i + 1
      else
         pphi(N_phi_t+2) = 0.0d0    !pphi5: C_L separable
      endif
      
      if (N_ssa_u>0) then
          pssa_62 = gr_Vpar(i)      !23
          pssa_65 = gr_Vpar(i+1)    !24
          pssa_65_t = gr_Vpar(i+2)  !25
          i = i + 3
      else
          pssa_62 = 0.0d0
          pssa_65 = 0.0d0
          pssa_65_t = 0.0d0
      endif      
      
      if (IHasHealthAge>=1) then
         pa_health1E = gr_Vpar(i)
         pa_health3B = gr_Vpar(i+1)
         pa_health3B_t = gr_Vpar(i+2)
         pa_health4D = gr_Vpar(i+3)
         pa_health4D_t = gr_Vpar(i+4)
         i = i + 5
      else
         pa_health1E = 0.0d0
         pa_health3B = 0.0d0
         pa_health3B_t = 0.0d0
         pa_health4D = 0.0d0
         pa_health4D_t = 0.0d0
      endif
      
      if (Ipt == 1) then  !part time option
         do j = 1, N_plpt 
            plpt_V(j) = gr_Vpar(i)
            i = i + 1   !index of next parameter if needed
         enddo  !do j
         do j = N_0, N_T
            plpt_t(j,1) = plpt_V(1)
            
            select case (Ipt_coef)
            case (1,21)   !polynomial, centered at N_0
                do ii = 2, 3
                    plpt_t(j,1) = plpt_t(j,1) + plpt_V(ii)*((dble(j-N_0)/10.0d0)**(ii-1))
                enddo !do ii
            case (2, 22)  !linear + alpha
                plpt_t(j,1) = plpt_t(j,1) + plpt_V(2)*(dble(j-N_0)/10.0d0)
                plpt_t(j,1) = plpt_t(j,1) + plpt_V(3)*pa0   !depends on pa0 as well
            case (3)  !centered at 45
                plpt_t(j,1) = plpt_t(j,1) + plpt_V(2)*(dble(j-45)/45.0d0) + plpt_V(3)*((dble(j-45)/45.0d0)**2)
            case (15)  !search center
                plpt_t(j,1) = plpt_t(j,1) + plpt_V(2)*((dble(j)-plpt_V(3))/plpt_V(3))  & 
                            & + plpt_V(4)*(((dble(j)-plpt_V(3))/plpt_V(3))**2)    
            case default
                if (j<=29) then
                    plpt_t(j,1) = plpt_t(j,1) + plpt_V(2)*(dble(29-j)/10.0d0)
                elseif (j>=50) then
                    plpt_t(j,1) = plpt_t(j,1) + plpt_V(3)*(dble(j-50)/10.0d0)
                else
                endif
            end select
            
            select case (Ipt_coef)
            case (21,22)
                plpt_t(j,2) = plpt_t(j,1) + plpt_V(N_plpt-1)
                plpt_t(j,3) = plpt_t(j,1) + plpt_V(N_plpt)
            case default           
                plpt_t(j,2:N_FS) = plpt_t(j,1)
            end select
            
            do ii = 1, N_FS
                plpt_t(j,ii) = 1.0d0/(1.0d0+exp(-plpt_t(j,ii)))
            enddo  !do ii
         enddo  !do j
      endif   !if (Ipt == 1) then  !part time option
      
      psigma_onjob = psigma  !psigma_onjob
      if (N_depreciation==2) then
          psigma_offjob = gr_Vpar(i)    !psigma_offjob (depreciation rate off job)
          i = i + 1
      else
          psigma_offjob = psigma
      endif
      
      if (gcI_robust == 22) then !22. depreciation is 0 for working period but still positive if not working
          psigma_onjob = 0.0d0
      elseif (gcI_robust == 23) then
          psigma_offjob = 0.0d0  !23. depreciation is 0 for not working but still positive if working
      endif
      
      if (N_pbequest == 2) then
          !setting Asset lower bound
          if (gcI_robust == 8 .OR. gcI_robust == 10 .OR. gcI_robust == 12) then  !exo
              rtemp2 = 25000.0d0 / grTotHours   !average annual income is 25k for high school graduates
          else
              rtemp2 = 500.0d0 / grTotHours
          endif
          sgridAnrlb%a(N_T_grid+1) = min(0.0d0, -pb2 + 10.0d0/grTotHours)    !grCF
          do j = N_T, N_0, -1
              ii = j - N_0 + 1   
              rtemp = max(0.0d0, dble(N_T-5-j)/dble(N_T-N_0)) * rtemp2      !gr_BL    !income in period git  
              sgridAnrlb%a(ii) = min(0.0d0, (sgridAnrlb%a(ii+1) - rtemp) * pr_inv)
              sgridAnrlb%a(ii) = sgridAnrlb%a(ii) * gr_LowerBound
              if (j >= N_0+1) sgridA%a(ii,1:VN_A(j)) = sgridA_BC%a(ii,1:VN_A(j)) + sgridAnrlb%a(ii)
          enddo !git
      endif  !if (N_pbequest == 2)
      
      if ((gcI_CF==88 .OR. gcI_CF==70) .AND. (new_id_mpi==0)) then
             write(*,'(A,I6,A,I2,A,F10.5,A,F10.5)') "iter#", gcI_itercount, &
                 & ', *** PARAMETERS in head slave of color ',color_mpi, ", pa0=", pa0, ", ppi=", ppi
      endif
      !============ parameters finished ================
      
      !==========some useful constants==========
      gr_atilde(1) = pa0          ! taste of leisure for singles
      gr_atilde(2) = pa0 + pa1    ! taste of leisure when spouse not working
      gr_atilde(3) = pa0 + pa2    ! taste of leisure when spouse working
      
      do i = N_0, N_T
         !taste shifter of consumption 
         gr_psi(i,1,0) = 0.0d0
         
         if (N_phi_t>0) then
             rtemp = dble(i-N_0)/1.0d1
             do j = 1, N_phi_t
                gr_psi(i,1,0) = gr_psi(i,1,0) + pphi(j) * (rtemp**j)
             enddo !do j
         endif
         
         gr_psi(i,2:3,0) = exp(gr_psi(i,1,0) + pphi(N_phi_t+1))   !married and leisure=0
         gr_psi(i,1,0) = exp(gr_psi(i,1,0))
         
         !for C_L nonseparable
         ! in separable case, pphi(N_phi_t+2) = 0.0d0 
         do j = 1, 3
             if (IC_L_nonseparable >= 1) then
                 gr_psi(i,j,1) = max(gr_psi(i,j,0)+pphi(N_phi_t+2), 1.0d-2)   !leisure=1
                 if (Ipt == 1) gr_psi(i,j,2) = max(gr_psi(i,j,0)+pphi(N_phi_t+2)*plpt_t(i,1), 1.0d-2)   !leisure=2 
             else
                 gr_psi(i,j,1) = max(gr_psi(i,j,0), 1.0d-2)   !leisure=1
                 if (Ipt == 1) gr_psi(i,j,2) = max(gr_psi(i,j,0), 1.0d-2)   !leisure=2
             endif
             gr_psi(i,j,0) = max(gr_psi(i,j,0), 1.0d-2)                   !leisure=0 
         enddo
         
         !Taste of leisure
         if ( (gcI_CF >= 46) .AND. (gcI_CF <= 50) ) then
            rtemp = dble(min(i,50)) !anticipated, taste for health does not increase with age after age 50
         else
            rtemp = dble(i)
         endif
             
         do j = 1, N_FS
            !health 1-Excellent 2-Good 3-Bad 4-Disabled 
            gr_leisure_taste(i,1,j) = gr_atilde(j) + (pa_health1E)        !health = Excellent (1)
            gr_leisure_condexp(i,1,j) = exp(gr_leisure_taste(i,1,j) + (pa_epsilon**2)/2.0d0)  !health = Excellent (1)
            
            gr_leisure_taste(i,2,j) = gr_atilde(j)                                            !health = good (2) BASELINE
            gr_leisure_condexp(i,2,j) = exp(gr_leisure_taste(i,2,j) + (pa_epsilon**2)/2.0d0)  !health = good  (2) BASELINE
            
            gr_leisure_taste(i,3,j) = gr_atilde(j) + (pa_health3B+pa_health3B_t*rtemp)        !health = bad (3)
            gr_leisure_condexp(i,3,j) = exp(gr_leisure_taste(i,3,j) + (pa_epsilon**2)/2.0d0)  !health = bad  (3)
            
            gr_leisure_taste(i,4,j) = gr_atilde(j) + (pa_health4D+pa_health4D_t*rtemp)        !health = disabled (4)
            gr_leisure_condexp(i,4,j) = exp(gr_leisure_taste(i,4,j) + (pa_epsilon**2)/2.0d0)  !health = disabled  (4)
         enddo !do j
      enddo !do i = N_0, N_T
      gr_leisure_condexp = gr_leisure_condexp * grUVmultiplier
      
      do j = 1, N_FS
         do i = 0, Ipt+1
            gr_bequest_AP(j,i) = 1.0d0 / (1.0d0 + (pbeta*pb1/gr_psi(N_T,j,i))**(1.0d0/petac))  !For the last period, with bequest
         enddo !do i
      enddo !do j

      gr_XI_sd = log(1.0d0 + pXI_sd**2) !variance of log(XI)
      gr_XI_mean = - gr_XI_sd/2.0d0   !mean of log(XI)
      gr_XI_sd = sqrt(gr_XI_sd) !sd of log(XI)
      
      gr_XI_xw(:,3)=exp(gr_XI_xw(:,1)*1.41421356237d0*gr_XI_sd+gr_XI_mean)
      
      git = N_T  !this is important. To set the backward induction at the beginning.
      git_pos = git - N_0 + 1
      
      return
  end subroutine sub_parse_parameters
  
  
!called in all nodes, once, in initialization
!A, AIME, H grids are same across het types.
subroutine sub_sgridA_set ()
  implicit none
  real(8) :: rtemp, Vrx2(2)
  integer(4) :: j1, j2
  character(100) :: filename
  

  !setting Asset lower bound
  if (gcI_robust == 8 .OR. gcI_robust == 10 .OR. gcI_robust == 12) then  !exo
      Vrx2(1) = 25000.0d0 / grTotHours   !average annual income is 25k for high school graduates
  else
      Vrx2(1) = 500.0d0 / grTotHours
  endif
  sgridAnrlb%a(N_T_grid+1) = min(0.0d0, -pb2 + 10.0d0/grTotHours)    !grCF
  do git = N_T, N_0, -1
      git_pos = git - N_0 + 1   
      rtemp = max(0.0d0, dble(N_T-5-git)/dble(N_T-N_0)) * Vrx2(1)      !gr_BL    !income in period git  
      sgridAnrlb%a(git_pos) = min(0.0d0, (sgridAnrlb%a(git_pos+1) - rtemp) * pr_inv)
  enddo !git
  sgridAnrlb%a = sgridAnrlb%a * gr_LowerBound
  
  do git = N_0, N_T
     git_pos = git - N_0 + 1  
     !set the asset grid, which is same for all (H,AIME,RECSS) at period git.
     if (VN_A(git) > 1) then 
         sgridcom%d1 = VN_A(git)
         allocate(sgridcom%a(sgridcom%d1))
         call ARRAY1D_INIT(sgridcom,VN_A(git),(/0.1d0,10.0d0/),1) !log linear
         rtemp = sgridcom%a(VN_A(git))-sgridcom%a(1)
     endif
     
     Vrx2(1) = sgridAnrlb%a(git_pos)  !asset natural lower bound
     
     if (Icol==0) then
        Vrx2(2) = 6.0d0 * gr_Amax%a(git_pos) !asset upper bound, min is 50  !130*10=1300 at age=49
     else
        Vrx2(2) = 9.0d0 * gr_Amax%a(git_pos) !asset upper bound, min is 50  !130*15=1950 at age=49 
     endif
    
     Vrx2(2) = Vrx2(2) + pinitA
         
     if (git>N_0) Vrx2(2) = max(Vrx2(2), 1.0175d0 * sgridA%a(git_pos-1,VN_A(git-1)))  !upper bound does not decrease
     
     if (gcI_robust==4) Vrx2(2) = Vrx2(2) * 0.5d0  !adjust the upper bound for gcI_robust==4 (pbeta = 0.95d0)
     
     if (git==N_0 .AND. VN_A(git)==2) then
         sgridA%a(git_pos,1) = Vrx2(1)
         sgridA%a(git_pos,VN_A(git)) = pinitA
     elseif (VN_A(git) > 2) then
         sgridA%a(git_pos,1:VN_A(git)) = Vrx2(1)+(Vrx2(2)-Vrx2(1))*(sgridcom%a(1:VN_A(git))-sgridcom%a(1))/rtemp
     else
         sgridA%a(git_pos,1:VN_A(git)) = pinitA
     endif
     
     if (VN_A(git) > 1) deallocate(sgridcom%a)
  enddo !git
       
  
  !==for debug only==
  if (0>1 .AND. gcI_itercount==0 .AND. new_id_mpi==1) then
     write(filename, "(A22,I1,A4)") "output/txt_MSM_Space_A",color_mpi,".txt"
     open(325,file=trim(cwd)//trim(filename)) 
     do j1 = N_0, N_T
        do j2 = 1, VN_A(j1)
           write(325, *) j1,",",j2,",",sgridA%a(j1-N_0+1,j2)
        enddo
     enddo   
     close(325)
     write(filename, "(A22,I1,A4)") "output/txt_MSM_Space_H",color_mpi,".txt"
     open(325,file=trim(cwd)//trim(filename)) 
     do j1 = N_0, N_T
        do j2 = 1, N_H
           write(325, *) j1,",",j2,",",sgridH%a(j1-N_0+1,j2)
        enddo
     enddo   
     close(325)
     write(filename, "(A25,I1,A4)") "output/txt_MSM_Space_AIME",color_mpi,".txt"
     open(325,file=trim(cwd)//trim(filename)) 
     do j1 = N_0, N_T
        do j2 = 1, VN_AIME(j1)
           write(325, *) j1,",",j2,",",sgridAIME%a(j1-N_0+1,j2)
        enddo
     enddo   
     close(325)
  endif  !if (1<0) then
  
  if ( (gcL_searchBC) .OR. (N_pbequest == 2) ) then
      sgridAnrlb_BC%d1 = N_T_grid +1
      allocate(sgridAnrlb_BC%a(sgridAnrlb_BC%d1))
      sgridAnrlb_BC%a = sgridAnrlb%a
      
      sgridA_BC%d1 = N_T_grid
      sgridA_BC%d2 = N_A
      allocate(sgridA_BC%a(sgridA_BC%d1, sgridA_BC%d2))
      sgridA_BC%a = sgridA%a
      do git = N_0, N_T
          git_pos = git - N_0 + 1
          sgridA_BC%a(git_pos,1:VN_A(git)) = sgridA_BC%a(git_pos,1:VN_A(git)) - sgridA_BC%a(git_pos,1)  !relative to the smallest element
      enddo !do git
  endif !if (gcL_searchBC)
  
end subroutine sub_sgridA_set

  
  subroutine normal_discrete (n_points, mu, sigma, x, w)
    !This is to discretize a standard normal distribution with mean mu and sd sigma
    !the grid is in x(:), and the weight (PMF) is in w(:)
    implicit none
    integer, intent(in) :: n_points
    real(8), intent(in) :: mu, sigma
    real(8), INTENT(OUT) :: x(:), w(:)
    integer(4) :: i
    real(8) :: test1, test2
    
    call Gauss_Hermite(x,w)
    x = x * 1.41421356237d0 * sigma + mu
    w = w / 1.77245385091d0    ! /sqrt(pi)
    

    !check mean and variance
    test1 = 0.0d0
    test2 = 0.0d0
    do i = 1, n_points
        test1 = test1 + w(i)*x(i)
        test2 = test2 + w(i)*x(i)*x(i)
    enddo
    
    return
  end subroutine normal_discrete

!======================================================================
!  utility
!======================================================================
  real(8) function ufunc_c(c_)
      !period utility function
      IMPLICIT NONE
      REAL(8), INTENT(IN) :: c_

      if (c_ <= 0.0d0) then
         ufunc_c = - NR_BIG2
      else
         if (abs(petac-1.0d0) < 1.0d-6) then
             ufunc_c = log(c_)
         else
             ufunc_c = ( c_ ** (1.0d0 - petac) ) / (1.0d0 - petac)
         endif
      endif

      ufunc_c = gr_psi(git, giFS, gI_p_L) * ufunc_c
      
      if (giRECSS==1 .AND. giRECSS_next==2 .AND. N_ssa_u>0) then !SSA at age git 
         if (git == 62) then
             ufunc_c = ufunc_c + pssa_62  
         else
             if (gcI_CF==9) then
                 if (git >= pNRA .AND. git <= 70) then
                     ufunc_c = ufunc_c + pssa_65 + pssa_65_t*dble(git-pNRA) 
                 endif
             else
                 if (git >= 65 .AND. git <= 70) then
                     ufunc_c = ufunc_c + pssa_65 + pssa_65_t*dble(git-65)
                 endif
             endif  !if (gcI_CF==9)
         endif      !if (git == 62)
      endif
      
      ufunc_c = ufunc_c * grUVmultiplier
      
      return
  END FUNCTION ufunc_c  

  !===========================================
  !  Map AIME to PIA
  real(8) function AIME2PIA( aime_ )
     IMPLICIT NONE
     real(8), INTENT(IN) :: aime_
     AIME2PIA = pPIAratio(1) * max(0.0d0, min(aime_, pPIAbp(1)))  & 
               & + pPIAratio(2) * min(pPIAbp(2)-pPIAbp(1), max(0.0d0, aime_ - pPIAbp(1))) &
               & + pPIAratio(3) * max(0.0d0, aime_ - pPIAbp(2))
  end function AIME2PIA
  !===========================================
  !===========================================
  !  Map PIA to AIME
  real(8) function PIA2AIME( pia_ )
     IMPLICIT NONE
     real(8), INTENT(IN) :: pia_
     PIA2AIME = max(0.0d0, min(pia_/pPIAratio(1), pPIAbp(1)) )  &
               & + min(pPIAbp(2)-pPIAbp(1), max(0.0d0, (pia_ - pPIAratio(1)*pPIAbp(1))/pPIAratio(2))) &
               & + max( 0.0d0, (pia_ - (pPIAratio(1)*pPIAbp(1) + pPIAratio(2)*(pPIAbp(2)-pPIAbp(1))))/pPIAratio(3) )
  end function PIA2AIME
  !===========================================
 
!======================================================================
!adjust the position of x in the state grids of the period t_
!ipos(1:3) pos of A/H/AIME in grid of sgridA/H/AIME
!======================================================================
  SUBROUTINE xPosadjust(t_, x, ipos)
    IMPLICIT NONE
    Integer, INTENT(IN) :: t_
    REAL(8), dimension(:), INTENT(INOUT) :: x  !(1A,2H,3AIME)
    Integer, dimension(:), INTENT(OUT) :: ipos
    integer :: t_pos

    t_pos = t_ - N_0 + 1
    ipos(1) = locate(sgridA%a(t_pos,1:VN_A(t_)), x(1))  !A
    if (ipos(1) == 0) then
       ipos(1) = 1
    elseif (ipos(1) == VN_A(t_) .AND. VN_A(t_)>1) then
       ipos(1) = VN_A(t_) - 1
       if (gcL_Iscapped) x(1) = sgridA%a(t_pos,VN_A(t_))
    endif
 
    ipos(2) = locate(sgridH%a(t_pos,:), x(2))  !H
    if (ipos(2) == 0) then
       ipos(2) = 1
    elseif (ipos(2) == N_H) then
       ipos(2) = N_H - 1
       if (gcL_Iscapped) x(2) = sgridH%a(t_pos,N_H)
    endif
    
    ipos(3) = locate(sgridAIME%a(t_pos,1:VN_AIME(t_)), x(3))  !AIME
    if (ipos(3) == 0) then
       ipos(3) = 1
    elseif (ipos(3) == VN_AIME(t_) .AND. VN_AIME(t_)>1) then
       ipos(3) = VN_AIME(t_) - 1
       if (gcL_Iscapped) x(3) = sgridAIME%a(t_pos,VN_AIME(t_))
    endif
  end subroutine xPosadjust
 
 
  !====================================================
  ! Bequest function
  !====================================================
  real(8) function Vbequest_A(A_)
  !the bequest value for the N+1 period, given asset A_
      IMPLICIT NONE
      REAL(8), INTENT(IN) :: A_
      if ( A_ <= -pb2 ) then
         Vbequest_A = - NR_BIG2 
      else
         Vbequest_A = pb1 * ( (pb2 + A_)**(1.0d0-petac) ) / (1.0d0 - petac)
      endif
      
      Vbequest_A = Vbequest_A * grUVmultiplier
  END function Vbequest_A  

  !====================================================
  ! tax code: output is the post-tax income
  ! gr_ftax_pbk(1) is 0.0d0
  !====================================================
  real(8) function ftaxcode(a_)
     implicit none
     real(8), INTENT(IN) :: a_
     integer(4) :: i

     do i=1, 6
        if ( a_ <= gr_ftax_pbk(git,i+1) ) then
           ftaxcode = gr_ftax_pbase(git,i) + gr_ftax_MPTR(git,i)*(a_-gr_ftax_pbk(git,i))
           exit
        endif
     enddo
     if ( a_ > gr_ftax_pbk(git,7) ) then
        ftaxcode = gr_ftax_pbase(git,7) + gr_ftax_MPTR(git,8)*(a_-gr_ftax_pbk(git,7))
     endif

     return
  end function ftaxcode
  
  !calculate the time elapsed
  !time_array(1) year
  !time_array(2) month of the year
  !time_array(3) day of the month
  !time_array(4) time offset with respect to UTC in minutes
  !time_array(5) hour of the day
  !time_array(6) minutes of the hour
  !time_array(7) seconds of the minute
  !time_array(8) milliseconds of the second
  function func_time_elapsed(dt0, dt1)
     implicit none
     integer(4), dimension(:), INTENT(IN) :: dt0, dt1
     integer(4), dimension(8) :: func_time_elapsed
     integer :: i
     
     func_time_elapsed(:) = dt1(:) - dt0(:)
     
     ! time_array(7) seconds of the minute, time_array(6) minutes of the hour
     do i = 7, 6, -1
        if (func_time_elapsed(i) < 0) then
           func_time_elapsed(i) = 60 + func_time_elapsed(i)
           func_time_elapsed(i-1) = func_time_elapsed(i-1) - 1
        endif
     enddo
     
     !time_array(5) hour of the day
     if (func_time_elapsed(5) < 0) then
        func_time_elapsed(5) = 24 + func_time_elapsed(5)
        func_time_elapsed(3) = func_time_elapsed(3) - 1
     endif
     
     !time_array(3) day of the month
     if (func_time_elapsed(3) < 0) then
        if (dt0(2)==2) then !Feburary
           !check if it is a leap year 
           if (mod(dt0(1),4)==0) then
              if (mod(dt0(1),100)==0) then
                 if (mod(dt0(1),400)==0) then    
                    i = 29 
                 else
                    i = 28 
                 endif
              else
                 i = 29 
              endif        
           else
              i = 28
           endif
        elseif (dt0(2)==4 .OR. dt0(2)==6 .OR. dt0(2)==9 .OR. dt0(2)==11) then
           i = 30
        else
           i = 31 
        endif  !if (dt0(2)==2)
        func_time_elapsed(3) = i + func_time_elapsed(3)
        func_time_elapsed(2) = func_time_elapsed(2) - 1
     endif
     
     !time_array(2) month of the year
     if (func_time_elapsed(2) < 0) then
        func_time_elapsed(2) = 12 + func_time_elapsed(2)
        func_time_elapsed(1) = func_time_elapsed(1) - 1
     endif
     
     !time_array(1) year
     
  end function func_time_elapsed
  
  subroutine create_directory( newDirPath )
    ! Author:  Jess Vriesema
    ! Date:    Spring 2011
    ! Purpose: Creates a directory at ./newDirPath

    implicit none

    character(len=*), intent(in) :: newDirPath
    character(len=256)           :: mkdirCmd
    logical                      :: dirExists

    ! Check if the directory exists first
    inquire( file=trim(newDirPath)//'/.', exist=dirExists )  ! Works with gfortran, but not ifort
 !   inquire( directory=newDirPath, exist=dirExists )         ! Works with ifort, but not gfortran


    if (dirExists) then
!      write (*,*) "Directory already exists: '"//trim(newDirPath)//"'"
    else
        mkdirCmd = 'mkdir -p '//trim(newDirPath)
        write(*,'(a)') "Creating new directory: '"//trim(mkdirCmd)//"'"
        call system( mkdirCmd )
        mkdirCmd = 'mkdir -p '//trim(newDirPath)//'/output'
        write(*,'(a)') "Creating new directory: '"//trim(mkdirCmd)//"'"
        call system( mkdirCmd )
    endif
  end subroutine create_directory

end module COMFUNCS
    